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ABSTRACT 

In this paper we carry out a quantitative analysis of the three-body systems and map 
them as a function of decaying time and intial conguration, look at this problem as an 
example of a simple deterministic system, and ask to what extent the orbits are really 
predictable. We have investigated the behavior of about 200 000 general Newtonian 
three body systems using the simplest initial conditions. Within our resolution these 
cover all the possible states where the objects are initially at rest and have no angular 
momentum. We have determined the decay time-scales of the triple systems and show 
that the distribution of this parameter is fractal in appearance. Some areas that appear 
stable on large scales exhibit very narrow strips of instability and the overall pattern, 
dominated by resonances, reminds us of a traditional Maasai warrior shield. Also an 
attempt is made to recover the original starting conguration of the three bodies by 
backward integration. We find there are instances where the evolution to the future 
and to the past lead to different orbits, in spite of time symmetric initial conditions. 
This implies that even in simple deterministic systems there exists an Arrow of Time. 

Key words: celestial mechanics - methods: A^-body simulations 



1 INTRODUCTION 

Dynamical interaction of multiple objects is common in 

■ physics and astronomy. It occurs over a huge range of scales 

■ from atoms in molecules to galaxies. If we know the in- 
[ teraction potential, then it is possible in principle to cal- 
culate the evolution of such a system. As is well known, 
systems of this kind tend to be chaotic, and so is the solu- 
tion of the three-bo dy problem (|Prigogine fc Stengerslll983 : 
lAarseth et"allll994l ). 

The two body problem in a Newtonian gravitational 
potential is simple enough to be solvable in a closed form. 
Adding a third body to the system complicates matters by 
prod ucing, in genera l, non-analytic solutions. Over a century 
ago (|Poincar j|l892l ') realized the chaotic nature of the few 
body system, but only with the advance of computer power 
and more sophisticated algorithms has it been possible to 
give more quantitative statements in the three bo dy problem 
jHeggielll975l : lMikkola|[l993 : lAarseth et al.lll994) . 

Here we introduce a simple way to illustrate and study 
the fractal nature of the general three body system. The 
method we describe provides an important tool for visual 
and qualitative studies of the three body problem. It allows 
us to handle and analyze all possible three body configura- 
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tions within our computational accuracy. The results visu- 
alised in this manner do not suffer of projection effects. With 
this method it is also easy to extend our study to include the 
influence of different initial conditions such as different ve- 
locities or masses. These extentions will help us understand 
the chaotic nature of the three body systems both in much 
greater detail and in a more comprehensive way. We also dis- 
cuss shortly the numerical sensitivity of general three-body 
problem and show that this may indicate the presence of the 
Arrow of Time within these systems. 



2 DECAY TIME-SCALE AND FRACTALITY 

We have considered the simplest initial state in an equal 
mass three-body system under Newtonian gravitation. The 
calculations start at rest and have no angular momentum. 
We have used in our calc ulations a fully chain-regu l arized 
Bulir sch-Stoer integrator (jMikkola fc AarsethI Il990l . Il993l , 
119961 ). 

To visualize our simulations we use an extended ver- 
sion of an homology map ( Agckian fc Anosoval Il967l ) . In 
this coordinate system one body resides initially at (0,0), 
the second one is at (1,0) and the third one is at (x,y). 
A shield shaped area contains all possible values of {x,y), 
spanning in abscissa from 0.0 to 1.0 and in ordinate from 
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Figure 1. The life time of gravitational three body systems containing all possible of orbits that were initially at rest. The initial 
configuration of the three body system is such that one star is located at (0, 0), the second one at (1, 0), and the third one at a coordinate 
{x,y). The colour pixel at {x,y) shows then the time elapsed from the initial configuration to the disintegration of the system under 
Newtonian gravitation. Red denotes short lived systems, decaying within one crossing time, orange systems decaying before two crossing 
times, etc. The dark blue colours represent systems with the longest lifespans (over 13 initial crossing times). The pixel size is 1.25- 10~^. 
The image of the complete shield is symmetric, because we start with all three objects at rest. Note the wide resonant bands causing a 
rapid disintegration of the system. Also note how narrower bands abound between wide bands. 



-73/2 to y/i/2 ^ 0.866. A three body confi euration can 
thus be unambiguously specified by a single point in this 
map. The unit time-scale of our (G — 1) system is t = 
^(i?/lAU)^/^ (M/Mq)"^''^ years, where R is the longest 
separation of the initial triangle and M is the mass of the 
objects. 

Except for a few dou ble or triple collision points 
l|Umehara fc Tanikawal |2000| ). all our three-body systems 
will eventually break up. The escaping body will have a non- 
negative total orbital energy. To define the breakup time of 
the system we also require that the escaping body must be 
at a distance of r ^ 1 from both remaining bodies. This 
is necessary because during strong interactions the total en- 



ergy of a body may, as expected, become momentarily highly 
positive. 

Because the escaping body often has only a mildly pos- 
itive energy and because the quantity we measure is time 
(the conjugate coordinate of energy in phase space), our 
measurement is in some respect analogous to a phase space 
cross-section of a triple system. 

The lifetime of the system may now be plotted as a 
function of the initial configuration. A complicated struc- 
ture emerges (Fig. [!}. As a curiosity, the overall pattern 
shows a nice r esemblance with a traditio nal African Maasai 
warrior shield llsee e.g. Huntingfordlll96ll ). Since we consider 
an equal mass system, there is a left-right symmetry in the 



Mapping the three-body system 3 




Figure 2. The life time of gravitational three body systems, a zoomed in image of Fig. [T] The pixel size here is 2.5 ■ 10~*. Colours as in 
Fig. [l] Note how this region appears quite 'stable' in Fig. [T] yet in this image it has broken into hundreds of narrow bands of instability. 
Note that the large yellow area breaks into three major parts. This is a universal feature for the loops. They correspond to the different 
object being ejected from the system. The point where these three parts meet is a triple collision point. The boundary of the yellow area 
is sharp from inside. If one approaches the boundary from the outside one encounters a number of orbits growing in density and stability 
as one approaches the sharp boundary. 



diagram. The non-zero initial velocity, eg. for the third body 
will in general break the up-down symmetry. 



2.1 Resonances 

On the largest scales the following regions of instability can 
be identified. A configuration initially close to an isosce- 
les triangle (1:1 resonance) is generally highly unstable (the 
vertical red band close to the vertical symmetry line of the 
map). Systems close to odd-digit orbital resonances (1:3, 1:5, 
1:7 etc) are highly unstable, because a very close approach 
takes place at the first close triple encounter. These are the 
bands growing denser towards the left and the right corners 
of the map. Note also the branching of the 1:7 and higher 
order resonances. Furthermore, resonances 1:2, 1:4 etc. are 
completely missing, because the phase in these cases is not 



favorable at the first close encounter to a rapid break-up of 
the system (jSaslaw fc Valtonenlll97j ). Systems breaking up 
due to the second close interaction show up as long orange 
arcs. The family of these arcs starts from the central part 
of the map and extends towards the edges of the map. A 
second system of arcs between the previous set and the 1:3 
resonance is also due to the odd resonances at the second 
encounter. A band due to similar resonances (nearly vertical 
bands) can be seen at the edge of the 1:1 resonance. 



All the major resonances appear to avoid a sea of appar- 
ent tranquility around the surroundings of {x, y) = (0.7, 0.5) 
and the corresponding points in the other quarters of the 
shield. 
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Figure 3. The life time of gravitational three body systems, a zoomed image of Fig. [2] showing the smaller details. The pixel size here is 
3.125- 10~^. The colour scale is linear and is extended to 20 crossing times. The dark blue colours show orbits that have not disintegrated 
in 20 crossing times. The light green band in the right upper corner of this image corresponds to a disintegration in about 7 crossing 
times. 



2.2 Local structure 

Locally the boundaries of these resonance regions are sharp 
from the inside. Approaching the boundary from outside 
one encounters ever denser loops and stripes. This structure 
appears selfsimilar and locally fractal. 

Areas of higher stability are also present in Fig [1] 
These fall into two categories. The first set are isolated deep 
trenches on large scales at about (0.70, 0.64) and (0.65, 0.15) 
and at respective symmetric points of the shield. The second 
set better visible on small scales are configurations just out- 
side some instability regions e.g. around the upper branch 
of the 1:7 resonance at (0.90, 0.23) or below the yellow eye 
at (0.75, 0.60) and the corresponding symmetric points. 

Figures [2] and |3] show progressive levels of zooming 
into these areas. Figure [3] represents a 5 by 5 pixel area 
in the lower left corner of Figure (2] The width of the strips 
seen in Figure [3] is comparable to the resolution element or 



3.125 • lO"''. These stripes are more easily seen if the im- 
age is viewed at a small angle. It is clear that a tiny change 
in the initial conditions may cause a substantial difference 
in the lifetime and in the evolution of the system. Further 
details of the orbital behavior within these areas is beyond 
the scope of this paper, except that we may state that these 
configurations appear to be more stable. 

2.3 Global fractality 

To investigate the ratio of stable (locally homogeneous) re- 
gions to the chaotic regions we compared, to, the ejection 
time of an orbit to tmcd, the median of 8 nearby orbits lo- 
cated on a rectangular grid with total width of 3.75 • 10"'^. 
The eight configurations are separated by e = ±1.25 • 10~^ 
in each of the x and y coordinates from the initial point. 

If l^mod — to\ ^ 0.5, we considered the point to have 
an locally homogeneous neighbourhood. Otherwise it was 
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Figure 4. The number fraction of recovered and unrecovered 
orbits as a function of required accuracy <5. At around 5 = 10~^ 
one come to the Umit of decimals in the result files. 




Figure 5. An example of a time symmetric orbit. The orbit is 
calculated first from free-fall position to escape (part not shown 
here) and then correctly back in time to the original free-fall situ- 
ation which is found at t = 3.6 and indicated by an arrow. Finally 
the integration was continued until the system breaks up again 
at t = 7.2. 



considered to have high sensitivity to initial conditions, a 
characteristic of chaos. Of the 21 699 cases compared in this 
manner 6212, or — 28.6 per cent were in locally homoge- 
neous. The remaining 71.4 per cent of positions turned out 
to be highly sensitive to initial conditions at this resolution 
implying that on the scale of e = 1.25 • 10~^ the fractal area 
covers about 71.4 per cent of all the systems. The impli- 
cation of this is that if one is able to determine the initial 
state of a triple system to an accuracy of e, one will, with 
the stated probability, be unable to tell the time-scale of the 
ejection of the third body to within half a crossing time. 

If our image is fractal one should find more and more 
(local) order as one decreases the resolution element. This 
will be difficult to do due to restrictions in computation 
time. We can increase the resolution and see if the amount 
of local homogeneity decreases in a self similar way. For bin 
sizes 1, 2, 3, 4, 5 and 10 times the original one we find 
a decrease in local homogeneity, pe, 28.6, 16.8, 14.9, 13.8, 
13.5 and 8.6 per cent, respectively. To within measurement 




Figure 6. An example of an orbit which is asymmetric with re- 
spect to time. This orbit is also traced back to its free-fall moment 
at t = 15.8 (marked with long arrow) but the motion after this 
differs from the previous history at t = 18.6 (marked with two 
arrows) leading in this case dramatically different break up. 



errors cx e~ . We find the same exponent if we study 
the detailed areas shown in Figures [2] and |3l The HausdorfT 
dimension ((Hausdorfiil919. ') of our image is thus 



1.5. 



3 TIME REVERSAL AND THE ARROW OF 
TIME 

The initial conditions we choose for our systems have the im- 
portant mathematical property that the orbits we calculate 
are time symmetric with respect to the origin. This implies 
that the triple systems we calculate are transient in the sense 
that the body expelled must have been once captured by the 
very same bodies forming the final binary system (and even 
in the very same phase but with an opposite sign!). Physi- 
cally this is not quite true, because even a small disturbance 
to the system would cause the system to loose this property. 
This also implies that a simple time reversal at the time of 
escape will not be a useful indicator of the errors involved 
in these calculations. 

It is possible that the initial conditions are not recov- 
ered by two reasons: the calculations may not be accurate 
enough, and errors accumulate to the extent where the ini- 
tial conditions are missed on the way back. The second pos- 
sibility is that the calculation is accurate enough, but the 
individual system is so sensitive that it is pra ctically impos- 
sible t o find the initial conditions again. See lAarseth et al.l 
(|l994h for a discussion of predictable and non-predictable 
orbits. 

To ensure independence from the method our calcu- 
lations were done with two different methods. The method 
uses a code consisting of a Bulirs ch-Stoer integrator with th e 
chain regularization algorithm (*Mik kola fc AarsethI Il993l ). 
The second code uses a simple leapfrog integrato r to log- 
arithmic Hamiltonian (|Mikkola fc Tanikawa|[l999l ). Results 
obtained are in agreement but in analysis we prefer the lat- 
ter code because of its constant stepsize. 

The orbits were calculated in few stages. First a for- 
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Figure 7. The homology map for most recovered orbits. The structure of the map is very similar to the one obtained by ejection times 
in Fig. [T] Here the colour indicates 5, the amount of time symmetry, for each orbit. Notice how this value tend to increase when initial 
configuration depart from the resonance bands. 



ward integration was carried out until one of the bodies es- 
capes. Here one should notice that criterion for recognizing 
an escape is not trivial. An escape is suspected if the ra- 
tio between the semi-major axis of one body orbit and the 
semi-major axis of the binary (i.e. two bodies close to each 
other) is large enough. We have selected this ratio to be 100. 
In addition to this requirement the system has to fulfil one 
of the two following conditions. An escape is identified if the 
third body is on an hyperbolic orbit and is moving away 
from the binary. Alternatively, if one of the bodies gets so 
far that the ratio exceeds ^ 2000 we consider this an escape, 
even though a proper hyperbo l ic orb it is not present. These 
escape criteria follow ^nosoval (|l986l ). which includes also a 
list of different possible escape criteria. We have also tried 
smaller values for the escape criteria with no apparent effect 
on our results. 

At the second stage the escape velocities were reversed 
and a reverse orbit was calculated until the system reached 
its initial conditions again. After that the simulation has 
carried on until the system breaks up again. 

The definition of when the initial conditions has been 
recovered depends on the accuracy of the computer and soft- 
ware. We find that by using double precision we can ex- 
pect results to be correct at least with accuracy ~ 10~^. To 
avoid rounding errors we consider only five significant digits. 
Therefore we require that the separation between original 



initial state and its recovered counterpart must be less than 
S = ±5.0 • 10~^. The choice of this limit affects naturally 
the ratio of the recovered and unrecovered orbits. The ratio 
is shown in Figure |4] where different values of 5 has been 
applied to the dataset. 

If one would have an infinite accuracy one would ob- 
tain an orbit where the motion is exactly the same before 
and after the free-fall moment - in the other words the or- 
bit is time symmetric. This is a fact which is based on the 
equations of motion. However, in most cases in our numeri- 
cal calculations this is not the case even though the free-fall 
moment can be recovered: most orbits actually forget their 
past after the free-fall moment producing an asymmetrical 
time evolution. This difference is illustrated with examples 
in FigureOfor time symmetric and Fig.|6]for time asymmet- 
ric orbits. We define a new quantity ^ as time ratio of the a 
recovered part to the orbit's time of disruption ttot- If ^ = 1 
the orbit is exactly the same in the past and in future from 
the free-fall moment. In Fig. [7] the value of ^ is presented 
for each orbit identified by its initial position on the homol- 
ogy map. The orbits marked by pink colour represent the 
totally symmetrical evolution while other colours indicate 
binned values of ^ for asymmetric orbits. In red pixels the 
length of the symmetric part is between 0.45 — 0.60 of ttot. 
In the yellow pixels symmetry covers 0.6 — 0.9 of ttot and 
at the blue pixels the orbit is over 90 but less than 100 per 
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Figure 8. The homology map for recovered orbits with ^ ^ 0.45. 
These orbits do not form any grouping as seen for higher values 
of e. in Fig. [3 



cent symmetric. The white areas - distributed quite equally 
everywhere on the map outside the resonances - loose their 
symmetry faster than within 0.45 of ttot or they represent 
cases where the orbit has not been recovered. The bins for 
^ in Fig. [7] are based on the distribution of orbits in Figure 
[9] which we will discuss later. 

According to Fig. [7| one will find the exactly time sym- 
metric behaviour within resonance areas where orbits have 
relatively short lifetimes. But is there any relation between 
the lifetime of the asymmetric orbits and their sensitivity to 
the numerical errors? This question is answered in Figure [9] 
where we present the amount of symmetry, f , as a function 
of disruption time of the system. It appears that the pop- 
ulation of asymmetric orbits are divided into two, partially 
overlapping, groups. Both of these groups have distinctive 
functional form ^i(*tot) = ijito"' + ^^'^ these sets appear 
to have quite sharp edges on the minimum side. Population 
I > 0.45) seems to be direct continuum of time symmetric 
orbits (C = 1) while the lower edge of the curve has approx- 
imately parameters ai = 0.5, ni = 0.7 and 6i = 0.47. These 
orbits appear to be located just outside the resonance ar- 
eas on the homology map and the value of ^ increases when 
the initial configuration moves farther from the resonance 
bands. 

Most orbits however belong to population II (0 < f < 1) 
which is quite equally distributed over the homology map, 
clearly further from the resonance bands. The lower limit 
of ^ in this population obeys approximate relation a2 — 
2.5, n2 = 0.9 and 62 ~ 0. On the homology map (Fig|8} this 
population is strongly intermixed with non-recovered orbits 
but the latter ones do not show any clear functional form 

fore = e(ttot). 

We have found out that even in quite a simple dynami- 
cal system such the classical equal-mass free-fall three-body 
problem, actually the Arrow of Time shows up. Considering 
the analytical point of view the autonomous form (i.e. no 
explicit time dependency) of the equations of motion and 
the lack of dissipation suggest that this arrow would not 
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Figure 9. The symmetricity, of the orbit as a function of the 
lifetime. The functional forms ^i(ttot) and 52(ttot) for two dis- 
tinctive populations of time asymmetric orbits are also shown 
(see text). Totally time symmetric orbits form short horizontal 
line at ^ = 1. 

be there if we could solve the problem exactly. However, in 
nature totally unperturbed systems do not exist and even 
a minute perturbation, just like a small computing error, is 
adequate to define a definite sense of time. 

How do we find out what is the forward sense of time 
in the three-body problem? The answer is provided by the 
concept of entropy again. Not the classical Boltzma nn en- 
tropy, but Kolmogorov- Sinai entropy (Kolmagorpj Il958l : 
ISinai|[l959l ). This concept of entropy applies to small sys- 
tems. It is related to the change of the volume of phase space 
occupied by systems which are initially very close to each 
other. It measures the loss of information about the initial 
conditions or the growth of disorder as time goes on. We 
have previously shown that Kolmogorov-Sinai entropy in- 
creases with time in the three-body proble m, and that it thus 
tells us what is future, and what is past l|Heinamaki et all 
119991 ). The Kolmogorov-Sinai entropy increases towards fu- 
ture time. 

There are other examples of determining the sense of 
time in the three-body problem. In three-body scattering a 
single body collides with a binary, and after some orbital 
evolution, one of the bodies escapes and leaves again a bi- 
nary. Let us say that the initial binary has zero eccentric- 
ity. The final binary h as an eccen t ricity e that follows the 
distibution /(e) = 2e ('jeansl ll919l : iHeggid [l975l ') . Thus the 
probability that the final binary escapes with zero eccentric- 
ity is zero. In a diagram dep icting such a scattering event 
(Isee e.g. Hut fc Bahcaiilll983l ) it is immediately clear which 
is the direction of orbital motion even if it is not indicated. 
Although some of our orbits perform time symmetric past 
and future they resemble a minority in our ensemble and 
it is somewhat reasonable to consider the Arrow of Time 
in the universe as a statistical quantity. One may also ask 
if these time symmetric orbits really appear at all in true 
nature where small perturbations are always present. 

The microscopic world consists of many small subsys- 
tems similar to the classical three-body system. It is our 
conjecture that our conclusions would apply also to such 
systems, and that Kolmogorov-Sinai entropy would actually 
provide the Arrow of Time for the microworld. 
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There have previously been a number of papers cosider- 
ing a relation between the Kolmogorov-Sinai entropy and 
the Boltzmann thermodynamic entropy. These studies, 
which have been done in other fields of physics. They 
have either tried to address the relation between these 
two quantities, albeit such r elation still remains u nknown, 
l|see e.g. Latora et al.l Il999l : iBar anger" et all l200lh or at- 
tempted to find a direct equivalence between the dy- 
namical instability a r id th e thermodynarn i c irreversibility 
jElskens fc Prigogind Il986l : lYukalovl I2OO3I : iRuiz fc Tsaliij 
I2OO7I ). If these relations indeed exist, they endorse our con- 
clusion. 



4 NOISE ESTIMATES 

Since noise estimates are difficult to perform in these kinds 
of systems. We followed the simpler practice of reducing the 
integration step by a factor of 2. We noticed only a slight 
sharpening in the boundaries, but no qualitative differences 
nor any unexpected quantitative differences. 

A second strong test we have performed is for the ro- 
bustness of the method. We ran a naive simulation with 
non-regularized equations of motion but with a short vari- 
able steps of integration. The code is very different in nature 
from our main code. Yet the results are very similar show- 
ing the main resonant structures. On small scales similar 
fractal-like patterns appear as with our main code. 



5 DISCUSSION 

The results from iHeinamaki et all (Il999l ) and lAnosoval 
l|l986t ). suggested that a homology triangle-like represen- 
tation is an ideal means to describe the evolution or the 
final solutions of a triple system, and that it would show 
structure both due to global and local sensitivity to initial 
conditions. The geometry of escape times which we have 
calculated shows rich structure featuring resonance bands 
at various levels and areas of apparent higher stability as 
well as areas with clearly a fractal structure. 

The time-scales in question for solar mass objects and 
an intial length scale of i? = 1 pc correspond to t ~ 
1.6 • 10^ years. The implication is that triple stellar systems 
forming initially from distant stars may still be in a multiple 
system after a few tens of million years (the blue areas in 
Figure [1] imply time-scales > 1.6 • 10* years). On the other 
hand one would expect stars being ejected by this method 
from young systems to form loose associations at the lat- 
est in about a ten million years. If the stars were closer 
[R = 0.1 pc) then the breakup times would be a factor 30 
shorter. On the other hand if the stars were not initially at 
rest then these time-scales are li kely to be somewhat longe r 
based on previous work (see e.g. ISaslaw fc ValtonenI |l97J); 
[Valtoncn & Karttuncn (200j)). 

Galaxies at distances of ~ 1 Mpc, and with masses M ~ 
50- IO^Mq have a typical time-scale of 10^'^'^ years, so these 
are expected to be still these initial orbits, even if they were 
near the resonance configuration. 
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